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We study the optimal packing of hard spheres in an infinitely long cylinder, using simulated 
annealing, and compare our results with the analogous problem of packing disks on the unrolled 
surface of a cylinder. The densest structures are described and tabulated in detail up to D/d = 2.873 
(ratio of cylinder and sphere diameters). This extends previous computations into the range of 
structures which include internal spheres that are not in contact with the cylinder. 



I. INTRODUCTION 

The dense packing of monodisperse (equal-sized) hard 
spheres in a cylinder has been found to produce a re- 
markable sequence of interesting structures as the ratio 
of the cylinder diameter to the sphere diameter is varied. 
We have explored computationally the densest of these 
packings in detail using simulated annealing, following 
the work of Pickett et aZ [1] , and the first four packings 
in this sequence are shown in Fig [I] We have reported 
extensive results, together with analytic approximations 
which serve to elucidate our findings [5]. The present 
paper and its future continuations amplify the previous 
reports by Mughal et al. [2] and by Chan (second author 
of the present paper) [3], and extend them in various 
directions. In particular, new results are provided for 
larger cylinder diameters and we include full details of 
all structures in a form suitable for future reference. 

We shall refer to these quasi ID packings as columnar 
crystals since they are periodic. Each structure can be 
assembled by stacking unit cells ad infinitum along the 
length of the cylinder with each subsequent unit cell ro- 
tated by the same twist angle with respect to the previous 
one. 

For smaller cylinder diameters, up to and beyond the 
range investigated by Pickett et al. [T] , the densest pack- 
ings only consist of spheres that are in contact with the 
surface of the cylinder. Eventually at larger cylinder di- 
ameters dense packings are found for which there are 
internal spheres that are not in contact with the surface. 
Nevertheless we extend the study of structures with only 
cylinder-touching spheres to larger diameters, for its own 
sake. As already explained in P], these structures are 
readily understood by recourse to a yet simpler problem, 
in which circular disks are placed on a cylinder. This can 
be fully worked out in analytical terms. 

In parallel with our computational study using sim- 
ulated annealing and a separate investigation by Chan 
[3] using sequential deposition, a number of experiments 
have been undertaken on the packing of solid spheres 
and small bubbles under gravity. Many of the simulated 
structures are observed. While we will have cause to 





FIG. 1. The first four densest sphere packings in a cylinder. 



mention these experiments, their presentation will be re- 
served for a further paper [4j. 

Some of the structures that we describe (particularly 
the more elementary ones) have been observed in bio- 
logical microstructures [S]. Columnar crystals have also 
been realised in numerous experimental contexts includ- 
ing foams (both dry [6Hl0] and wet [4]) in tubes, colloids 
in micro channels f llH15| . and fuUerenes in nanotubes 



We have not attempted to rigorously prove that the 
densest packings identified by simulations are indeed the 
densest possible. This ought to be achievable for some of 
the simplest ones. 



II. MOTIVATION AND BACKGROUND 

The identification of dense packings has always played 
an important part in condensed matter physics and phys- 
ical chemistry [T^. The entities which are to be packed 
are often spherical, or well approximated as such, and 
they may be monodisperse. In the idealised case of 
hard spheres (with no interaction when not in contact) 
the problem of finding the maximum density with given 
boundary conditions constitutes a classic mathematical 
challenge. If posed for an unbounded packing, it is long 
associated with the name of Kepler. It is still debated, 
in terms of full mathematical rigour, although there is 
no room for practical doubt about the nature of densest 
packings (fee etc) [20] . 

Here we are concerned with the case in which spheres of 
diameter d are to be contained in an unbounded cylinder 
of diameter D. We shall pursue it up to D/d = 2.873 
which represents the limit of the capability of our present 
computational resources and methods. 

Our primary objective is to identify the structure with 
the largest value of <&, defined as the fraction of volume 
occupied by the spheres. 

There are close connections between this topic and that 
oi phyllotaxis, a subject arising out of biology and having 
to do with the dense arrangement of similar units on the 
surface of a cylinder, exemplified by pine cones, pineap- 
ples or corn cobs [H]. This connection will be explored 
in detail here. 



III. SOME SIMPLE COLUMNAR CRYSTALS 
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For certain specific values of D j d the optimal packing 
structure can be easily surmised. These structures are di- 
rectly related to the analogous two-dimensional problem 
of finding the smallest diameter circle into which N non- 
overlapping circles, each of diameter d, can be packed. A 
description of these special cases will serve to illustrate 
the general problem that we will address. 

The first five of these special solutions, labelled Cat, 
are shown schematically in Fig [2] and (for A^ > 2) consist 
of disks placed at the vertices of a regular polygon. The 
diameter of the enclosing circle is given by. 
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where D''{2) = 2d, 
D%5) = 2.7013d. 

Replacing the circles with TV spheres we define a unit 
cell which can be repeated along the length of the tube as 
follows. Each successive layer, indicated by alternating 
red and green spheres, can be generated by translating 
the previous layer through a distance 
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FIG. 2. Top: The first five solutions of the circle packing 
problem. Bottom: Five columnar structures corresponding 
to the first five solutions of the circle packing problem. The 
unit cells are shown in red/green. 



along the tube and rotating it by a twist angle a''{N) — 
tt/N. Using the above formula the volume fraction <I> 
of these simple columnar crystals can easily be com- 
puted, $(1) = 2/3, $(2) = 0.4714, $(3) = 0.5276, 
$(4) = 0.5441, $(5) = 0.5370. 

We shall label these simple columnar crystals with the 
notation Cn to indicate that their structure can be de- 
rived in a simple way from the circle packing problem. 



IV. SIMULATION 

For general values of D/d, the optimal packing struc- 
ture cannot be guessed easily and we must turn to heuris- 
tic methods. Our search method is confined to structures 
that are periodic in the following sense. There is a primi- 
tive cell, of length L, containing N spheres, the structure 
being generated from this by the screw operation of (i) 
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FIG. 3. Volume fraction of the densest packing as a function of D/d. Discontinuities in the derivative are indicated by vertical 
lines. Red dotted lines indicate solutions to the circle packing problem. The main graph shows results in the range 1.8 < D/d < 
2.71486, the inset are continuations of the same graph showing the regions 1 < D/d < 1.8 (left) and 2.71486 < D/d < 2.873 
(right). Structures above D > 2.71486 include internal spheres while those below do not, the division between these two regions 
is denoted by a heavy blue dashed line. 



translation along the cylinder axis by nL (where n is any 
integer) combined with (ii) rotation about the axis by an 
angle na. This screw operation represents the underly- 
ing symmetry of columnar crystals (which are not to be 
confused with columnar phases, originating in the study 
of liquid crystals and related instead to the packing of 
columns). 

Our primary method of simulation is based on the well 
known approach of simulated annealing. This provides a 
reasonably exhaustive and unbiased search for maximum 
density. Appendix A summarises the technical details of 
the present application (such as annealing schedules). The 
search procedure described in Appendix A looks for the 
minimum possible value of L, for a given N, treating 
a and the sphere positions as variables. It does this for 
increasing values of N until the (likely) optimal structure 
is apparent, or further computation is not practical. 

In the previous simulations of Pickett et al. [T] , a pe- 
riodic boundary condition without any twist was used. 



If for example the densest structure has a twist angle a 
equal to 7r/5, then ten cells would combine to satisfy the 
simpler boundary condition. 



V. NUMERICAL RESULTS 

In this section we present the full range of results for 
several significant properties as a function of D/d. These 
are the volume fraction and chirality, for the densest 
structures that we have found. The following figures and 
Table U summarise the main results. 

Table H] enumerates all the densest packings observed, 
from our numerical work, in the range 1 < D/d < 2.873. 
The table lists the range of D/d over which each struc- 
ture is observed, the number of spheres in the primitive 
unit cell (from which the extended structure can be con- 
structed) and the average number of contacts per sphere. 
The table also classifies the structures into two groups: 
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FIG. 4. The numerically computed derivative (finite differ- 
ence) of the volume fraction curve, as shown in Figls] 




FIG. 5. A schematic representation illustrating variations in 
the volume fraction of the densest packing arrangements as a 
function of D/d (see text) 



those which arc chiral and those which are not (achiral) 
- see below for a more detailed discussion of chirality. 

The break in Table [ij between structures 32 and 33, 
highlights the fact that beyond D/d = 2.7379 the densest 
structures are those which contain spheres not in contact 
with the cylindrical boundary. The unit ceU for packings 
with internal spheres are described using the notation 
(771, n), where the first integer is the number of internal 
spheres and the second is the number of spheres which are 
in contact with the cyhndrical boundary. So for example 
in the case of structure 30 the unit cell contains 6 spheres 
of which 5 are in contact with the cylindrical boundary 
and 1 is not. 



explored. By simple finite difference we approximate the 
derivative of the volume fraction as a function of D/d, 
as shown in Fig|4j This is to clearly identify the singular 
behaviour discussed below. 

The vertical lines indicate a discontinuity in the deriva- 
tive of the volume fraction; there is either (i) only a 
sudden changes in the derivative, or (ii) the simulta- 
neous presence of such a sudden change and a square- 
root singularity in the derivative, the two cases of which 
are denoted by vertical blue-dashed and red-dotted lines, 
respectively. Note the square-root singularities in the 
derivative coincide with the C'n (circle packing) struc- 
tures. 

Let us discuss this behaviour by reference to Fig (Is]) 
which is a simplified cartoon of the observed changes in 
the volume fraction as a function of D/d. Qualitatively, 
the variation of the maximum density takes the following 
form: its dependence on D/d is everywhere continuous 
while the derivative is piecewise continuous. We distin- 
guish two types of singular points, as in Fig (pi), at which 
the derivative changes. 

The first occur when a maximum number of contacts 
is reached. These points correspond to highly symmetric 
structures, such as the C„ structures described above. At 
such points - represented by the vertices labelled a and 
c, in the inset - the previous trend of structural change 
cannot be continued; in effect, the structure has jammed 
due to the formation of new contacts. In other words: 
if the separations of contacting sphere centres is to be 
maintained, the system is overconstrained and no such 
solution exists beyond this point. Instead, some existing 
contacts are released and a new structural trend pro- 
ceeds. The structure itself varies continuously through 
the singular point, with a downward change in the deriva- 
tive of the density, corresponding to the line segments ab 
and cd in the cartoon. 

At the second kind of singular point the structure 
is simply overtaken by another more dense packing. 
Here the structure itself changes discontinuously, and the 
derivative obviously must change in a positive sense. In 
the inset this transition corresponds to the points labelled 
b and d where, for example, the line segment be illus- 
trates the increasing trend in density immediately fol- 
lowing such a transition. The dashed lines indicate the 
continuation of the previous structure which now has a 
lower density compared with the optimal packing. 



Volume fraction 



Fig [3] presents the maximum density found by our pre- 
scribed procedure, for the full range of D/d that we have 



These remarks are based on observation of the results. 
We have for example no proof that the change of deriva- 
tive at points of the first type is always negative, al- 
though we can offer a proof that the dependence of (max- 
imum) density on D/d cannot have any discontinuities at 
which the density undergoes a finite upward or downward 
change. This is discussed in detail in appendix C. 



Structure 


Range 


Number of 
spheres in 
Unit Cell 


Notation 


Average 
Contact 
Number 


Description 


chirahty 


l(Ci) 


D/d^ 1 


1 


straight chain 


2 


- 


achiral 


2 


1 < D/d< 1.866 


1 


zigzag 


2 


- 


achiral 


3 


1.866 < D/d < 1.995 


1 


twisted zigzag 


4 


- 


chiral 


4 


1.995 < D/d < 2.0 


2 


(2,1,1) 


4 


hne-slip 


chiral 


5(C2) 


D/d = 2.0 


2 


(2,2,0) 


5 


maximal contact 


achiral 


6 


2.0 < D/d < 2.039 


2 


(2,2,0) 


5 


hnc-shp 


chiral 


7 


D/d = 2.039 


1 


(3,2,1) 


6 


maximal contact 


chiral 


8 


2.039 < D/d < 2.1413 


2 


(3,2,1) 


5 


hne-shp 


chiral 


9 


2.1413 < D/d < 2.1545 


3 


(3,2,1) 


16/3 


hne-shp 


chiral 


10 (Cs) 


D/d = 2.1547 


3 


(3,3,0) 


6 


maximal contact 


achiral 


11 


2.1547 < D/d < 2.1949 


3 


(3,3,0) 


16/3 


hne-shp 


chiral 


12 


2.1949 < D/d < 2.2247 


2 


(3,2,1) 


5 


hne-shp 


chiral 


13 


D/d = 2.2247 


2 


(4,2,2) 


6 


maximal contact 


achiral 


14 


2.2247 < D/d < 2.2655 


2 


(4,2,2)\(4,2,2) 


5 


line-slip 


chiral 


15 


2.2655 < D/d < 2.2905 


3 


(3,3,0) 


16/3 


hne-slip 


chiral 


16 


D/d = 2.2905 


1 


(4,3,1) 


6 


maximal contact 


chiral 


17 


2.2905 < D/d < 2.3804 


3 


(4,3,1) 


16/3 


hne-slip 


chiral 


18 


2.3804 < D/d < 2.413 


4 


(4,3,1) 


11/2 


hne-slip 


chiral 


19 (C4) 


D/d = 2.4142 


4 


(4,4,0) 


6 


maximal contact 


achiral 


20 


2.4142 < D/d < 2.4626 


4 


(4,4,0) 


11/2 


hne-slip 


chiral 


21 


2.4626 < D/d < 2.4863 


3 


(4,3,1) 


16/3 


hne slip 


chiral 


22 


D = 2.4863 


1 


(5,3,2) 


6 


maximal contact 


chiral 


23 


2.4863 < D/d < 2.5443 


3 


(5,3,2) 


16/3 


hne-slip 


chiral 


24 


2.5443 < D/d < 2.5712 


4 


(4,4,0) 


11/2 


hne-slip 


chiral 


25 


D = 2.5712 


1 


(5,4,1) 


6 


maximal contact 


chiral 


26 


2.5712 < D/d < 2.655 


4 


(5,4,1) 


11/2 


line-slip 


chiral 


27 


2.655 < D/d < 2.7013 


5 


(5,4,1) 


28/5 


hne-slip 


chiral 


28 (Cs) 


D/d = 2.7013 


5 


(5,5,0) 


6 


maximal contact 


achiral 


29 


2.7013 < D/d < 2.71486 


5 


(5,5,0) 


28/5 


hne-slip 


chiral 


30 


2.71486 < D/d < 2.7306 


6 (1,5) 




26/6 


-(1,5) 


chiral 


31 


D/d = 2.7306 


6 (1,5) 




40/6 


maximal contact 


chiral 


32 


2.7306 < D/d < 2.74804 


6 (1,5) 




26/6 


+(1,5) 


chiral 


33 


2.74804 < D/d < 2.8211 


11 (1,10) 




40/11 


-(1,10) 


see Sect. VIII 


34 


D/d = 2.8211 


11 (1,10) 




60/11 


maximal contact 


achiral 


35 


2.8211 < D/d < 2.8481 


11 (1,10) 




58/11 


+ (1,10) 


chiral 


36 


2.8481 < D/d < 2.8615 


7 (2,5) 




32/7 


-(2,5) 


chiral 


37 


D/d = 2.8615 


7 (2,5) 




34/7 


maximal contact 


chiral 


38 


2.8615 < D/d < 2.8711 


7 (2,5) 




30/7 


+ (2,5) 


chiral 


39 


2.8711 < D/d < 2.873 


15 (2,13) 




72/15 


-(2,13) 


chiral 


40 


D/d = 2.873 


15 (2,13) 




90/15 


maximal contact 


chiral 



TABLE I. Specification of densest structures, up to D/d = 2.873. Bold numerals designate the line-slip type (see text section 
VILA). The break in the table denotes the point beyond which packings with internal spheres are observed. 



B. Chirality 

A chiral object is one which cannot be superimposed on 
its mirror image (or inverse)by translations and (proper) 



rotations. The question of chirality is interesting in the 
present context, for example in relation to designing chi- 
ral molecular filters that will discriminate between enan- 
tiomers. 



X 0.2 




FIG. 6. A plot of the chiral index x(D/d), which measures the degree of mismatch between a packing and its mirror image 
(see text). Vertical lines indicate structural transitions as in Fig 2. The main graph, on the left shows results in the range 
1.8 < D /d < 2.71486 while the graph on the right is a continuation but focuses on the region 2.71486 < D/d < 2.873. 



According to this definition an object is either chiral 
or it is not (it is achiral). The reported structures are 
classified accordingly in Table |l] 

But chirality may be manifested in physical properties 
to a greater or lesser extent: the degree of chirality of 
the structure itself is a tempting concept. While it may 
be useful to think of this in practice, there is no unique 
quantitative definition that will have general relevance. 
Any appropriate property could be used as an index of 
chirality. Nevertheless one may offer a working definition 
for use in relation to simulations, as Pickett et al did 1]. 
Here we will employ another definition that is perhaps 
more transparent. 

We define a chirality index x in terms of the degree 
to which an object can be superimposed upon its mirror 
image. For a given columnar packing, our method is to 
start with the sphere centres and generate a mirror image 
of the structure by reflecting the coordinates of the cen- 
tres in the x-y plane, i.e. the cross-sectional plane of the 
cylinder. For each sphere in the packing we compute its 
distance to the nearest sphere in the mirror image. The 
overlap function is defined as the sum of these distances 
divided by the number of spheres in the original packing. 
Clearly when a packing can be completely superimposed 
upon its mirror image the overlap function vanishes. 

The computational challenge then is to find the ar- 
rangement of the mirror image (by rotation and transla- 
tion) that minimises the overlap with respect to the orig- 
inal structure. This is done in a straightforward manner 
by simulated annealing. Clearly for chiral structures the 
overlap cannot vanish and by plotting the overlap func- 
tion we have a measure x of tho chirality as a function 
oi D/d, as in Fig([6]). 



VI. DETAILS OF STRUCTURES FOR D/d UP 
TO 2.71486 

Up to D/d = 2.71486, all the structures that are found 
are of a special character, helpful for their interpreta- 
tion. Every sphere is in contact with the cylinder sur- 
face. Therefore all of their centres lie on an inner cylin- 
der, of diameter D-d, and they touch another cylinder of 
diameter D-2d. They may therefore be considered as the 
densest packings of spheres on a cylinder. 

Clearly inner spheres must appear in the optimal pack- 
ing when D-2d is greater than d, that is D > 3.0, in fact 
they are first found at D/d > 2.71486, and we reserve 
them for a separate discussion in the next section. We 
also treat structures for low D/d separately in section 

rvTci 



A. MELximal contact structures 

We have already noted the special points at which 
the structure changes continuously, with the formation of 
new contacts and the breaking of old ones. We have pre- 
viously called them symmetric structures, but this term 
is really only meaningful in the disk-packing analysis be- 
low. Here we shall apply the term maximal contact in- 
stead. For D/d > 2, all of these structures are depicted 
in Fig (It]) and labelled using the phyllotactic notation 
explained in Appendix B. 

They include the simple C„ packings which we have al- 
ready described and are labelled n, n, in the phyllotactic 
notation, see Appendix B. This relates to the pattern of 
sphere centres and contacts, "rolled out" on to a plane. 



(2,2,0) 



(3,2,1) 



r (3,3,0) 




FIG. 7. Maximal contact structures, with the corresponding "rolled-out" pattern of contacts. The vector V 
(corresponding to the perimeter of a cylinder cross-section) is indicated by the red arrows. 



and identified with rhombic or triangular lattices. This 
does not apply to the first of the structures, the straight 
chain, and is rather confusing if applied to the next two, 
hence we use it only for the subsequent structures. 

We begin with the case of D/d ~ 2: this is the close- 
packed, achiral, structure C2 previously described. Each 
sphere has five contacts (not including contact between a 
given sphere and the cylinder); these consist of a contact 
with the neighbouring sphere in the same unit cell and 
four contacts with spheres in adjacent cells. All of the 
rest of the maximal contact structures so far considered 
are composed of spheres with six contacts. Packings of 
the type C„ are achiral while the rest are chiral. 



along one line separating two of the spiral chains of the 
symmetric structure, and a consequent relative slip of the 
two sides. Such line slip arrangements are a feature of 
hard sphere packings and are evident in the rolled out 
patterns of Fig [8) 

There are in general three possible choices for the di- 
rection of the line slip and the results of Fig [8] correspond 
to those that maintain the highest density as in Fig [3] 

We continue to higher values of D/d in section VIIL 



C. Structures below D/d = 2.0, and the square root 
singularity 



B. Line slip structures 

We now turn to the intervening structures, which are 
modified versions of the maximal contact packings, ad- 
justed to fit around the cylinder. Such staggered helices 
were first observed by Pickett et al in their simulations. 
We label these line slip structures to indicate that the 
modification consists of the release of half of the contacts 



We return to examine the structures for low D/d, 
which are somewhat different. The reason for such a 
difference will become clearer when we explore the rela- 
tionship to disk packing, in a later section. These first 
four structures are shown in FigjT] 

Structure 1 is the elementary case of a straight chain of 
spheres, like peas in a pod, as shown in Fig[l] This trivial, 
achiral, structure has a volume fraction of $(1) — 2/3. 
Similarly structure 2 is obvious and consists of a zigzag 
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FIG. 8. Line slip structures at arbitrarily chosen values of D/d within their respective ranges, with corresponding rolled out 
patterns of contacts. The vector V is indicated by an arrow. 



planar arrangement of spheres such that the change of 
the azimuthal angle from one ball to the next is equal to 
tt; each sphere makes contact with one sphere from above 
and one from below. When the zigzag packing, i.e struc- 
ture 2, encounters additional contacts between its second 
neighbours at D/d = 1.886 (so that each sphere now has 
four contacts), it is forced to take the form of an increas- 
ingly twisted spiral (structure 3). By direct numerical 
calculation one may follow this to its end point at struc- 
ture 7, but structure 6 intervenes with higher density, so 
that there must be a transition below D/d = 2.0. This 
is found at D/d=1.990, strikingly close to 2.0, beyond 



which the densest packing is a line-slip modification of 
structure 6. The smallness of the interval in which this 
is found is largely due to the existence of a square-root 
singularity. 

The variation of volume fraction, as any of the C„ or 
n, n, structures is approached with increasing D/d, ex- 
hibits a square-root form, whereas it is linear in other 
cases. That is 



<i)o - •& 



Do 
do 



1/2 



(3) 



where the quantities with a subscript belong to a C; 



9 



or n, n, structure. The derivative plot shown in Fig 2 is 
intended to make this singular behaviour more evident. 
The square-root arises from the special symmetry of these 
achiral structures. 



VII. RELATIONSHIP TO CIRCULAR DISK 
PACKINGS ON A CYLINDER 

Many aspects of these results are made more under- 
standable by recourse to the packing of circular disks on 
a cylinder. We have already noted that the same se- 
quence of symmetric structures is found. The similarity 
extends further to the line-slip structures and some of 
the details of the analytic form of the curve in Fig (Is]), 
not yet noted. 

This qualitative correspondence attracts our attention 
to the cylindrical disk packing problem, which we pursue 
below. 



A. Disk Packings 

Densest packing of circular disks in a plane places their 
centres on a triangular lattice where each disk has six 
contacts; this was rigorously proved about a century ago 
[l9] . Cylindrical surface packing of disks with the same 
density is generated by rolling this pattern on to a cylin- 
der, when possible. 

To see how this can be achieved seamlessly let us define 
a periodicity vector V between any pair of disk centres of 
the triangular lattice packing as shown in Fig ([9]) a. The 
region between the start and end of V, which is bounded 
by lines perpendicular to V, can then be excised and 
wrapped onto a cylinder whose circumference is equal 
to the length of V. The resulting structure is a dense, 
homogenous packing - in the sense that all disk sites are 
equivalent - which we call a symmetric packing. Any such 
symmetric packing is characterised by this periodicity 
vector. In the phyllotactic scheme V can be defined by 
a set of three ordered positive integers {I — m + n,m,n), 
as discussed in appendix B (a simple operational method 
to assign indices is to count the number of lattice rows 
that cross the periodicity vector in the three directions, 
as shown in Fig (|9])a and Fig ^c). 

For other values of the cylinder circumference, the sym- 
metric packing may be distorted in some way to wrap 
around the cylinder. The most obvious adjustment is an 
affine transformation of the triangular lattice, as shown 
in Fig Mb. We illustrate the effect of the transforma- 
tion by reference to the shaded triangle in Fig Ma. The 
affine transformation distorts the equilateral triangle into 
an isosceles triangle by keeping the length of two adjacent 
sides fixed and varying the third; consequently, the disk 
centres form a rhombic lattice which has a lower den- 
sity (area fraction) compared with the triangular lattice 
since each disk now has only four contacts (the lowest 
area fraction corresponds to a square lattice). The trans- 



formation may be adjusted to make the length of V equal 
to ttD. The resulting pattern can be wrapped onto the 
surface of the cylinder. We may call these asymmetric or 
affine lattice packings. 

If D/d is varied the disk centres of such a structure are 
eventually brought back into coincidence with the sites 
of a triangular lattice, see Fig Qc. In this manner the 
strained structure proceeds from one symmetric packing 
to another. Since there are three possible choices for 
the affine transformation, the rules for this process - as 
reported previously [2] - are, when applied to the second 
and third phyllotactic indices, as follows, 

(to, n) — > (m — n, n) 
[m, n) — > [m + n, m) 
{m, n) — > (m + n, n), 

where the new phyllotactic indices may have to be rear- 
ranged into descending order. Intermediate asymmetric 
packings may be labelled using rhombic notation [p, q] , 
where the ordered indices p > q are the indices common 
to both the initial and final symmetric states. 

We illustrate the use of these rules with an example. 
In the case of the symmetric packing [5, 4, 1] (illustrated 
in Fig (|9])a) the above rules yield. 



(4,1) 


-^ 


(3,1) 


(4,1) 


-^ 


(5,4) 


(4,1) 


-^ 


(5,1) 



Thus the symmetric packing [5, 4, 1] is connected to 
[4, 3, 1] via the rhombic structure [4, 1], to [9, 5, 4] via [5, 4] 
- which are shown in Fig (l9])c and Fig (l9])b, respectively, 
and to [6, 5, 1] via [5, 1]. 

Although such asymmetric packings are the simplest 
type of dense structure, intermediate between two sym- 
metric packings, they are not the densest. As reported 
previously [5] , asymmetric packings are in general super- 
seded by another type of packing involving an inhomo- 
geneous shear of the symmetric lattice in which there is 
a localised strain or slip along a line (and its periodic 
replicas) . 

We can describe such line-slip packings with reference 
to the symmetric packing [5,4, 1]. Again there are three 
possibilities and we describe each of these in turn. In 
the first case, four lattice rows cross the periodicity vec- 
tor in the horizontal direction, as shown in Fig Qa. We 
allow the final, fourth, row to "slide-over" the previous 
row until the disk centres are once again arranged into a 
triangular lattice. Thus as illustrated in Fig ^, the line 
slip in question is intermediate between the the symmet- 
ric packing [5,4, 1] and [6,4,2]. By allowing the disks in 
the final row to slide in the opposite manner we find that 
the symmetric packing [5,4,1] is connected to [4,4,0]. 
In the second case five lattice rows cross the periodic- 
ity vector. Four rows are held fixed while the final layer 
of disks slides over the penultimate row - thus [5, 4, 1] is 
connected to [5,3,2] and [5,5,0]. In the third case only 
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FIG. 9. (a) symmetric packing [5,4,1], the black arrow represents the periodicity vector; (b) the square packing [5,4]; (c) the 
symmetric packing [9,5,4] which is connected by an affine shear to [5,4,1]; (d & e) show a hne slip connecting [5,4,1] to the 
symmetric packing [6,4,2] (f). 



a single lattice row crosses the periodicity vector and we 
find that [5, 4, 1] is connected to [6, 5, 1] and [4, 3, 1] 

This localised shear allows the length of V to vary 
continuously. The disks involved in the line-slip have five 
contacts while the rest of the structure remains symmet- 
ric and close-packed, with six contacts for every disk. 
Clearly the area fraction has a maximum, corresponding 
to the symmetric packings, while the intermediate line- 
slip packings have a lower value. As reported previously 
[2], there is a simple rule for the close-packed structures 
that are the end points of line-slip solutions 

(1) {m,n) -^ (m + l,n) or (rn — l,n) 

(2) (tti, n) — > (?Ti, n -f 1) or (rn, n — 1) 

(3) (m, n) — ^ (?7i + 1, n — 1) or (m— l,n-f-l), 

where the leading numbers denote the direction ai, a2 
or §3 of the line-slip. Again the above rules apply to 
the second and third phyllotactic indices of a given close 
packed structure and keep either n, m or I constant. For 
example, using the above rules, the symmetric packing 
[5, 3, 2] is connected by a line slip along ai to [6, 4, 2] or 
[4,2,2], a line slip along a2 yields [4,3,1] or [6,3,3] and 



a hne slip along a^ yields [5, 4, 1] or [5, 3, 2] (note in the 
third case - along as - a rearrangement of the phyllotactic 
indices into descending order was necessary and the new 
structure is the same as the initial structure). 

In TableUwe use bold numerals to denote the direction 
of the line slip. In the table we only have cause to mention 
the first column of line-slip structures, that is: we denote 



(1) [l,m,n] 

(2) [l,m,n] 

(3) [l,m,n] 



^ -|- 1, TO + 1, n] by [Z,77i, n] 
^ + l,TO,,n + l] by [Z,m, n] 
I, m + 1,71—1] by [l,?7i,n], 



A full derivation of such connection rules for both the 
asymmetric and line-slip packings will be given in a fu- 
ture publication. 

Fig (10 1 presents analytical results for the density of 
these disk packings. Here we see for the first time some 
of the structures of lower density, as well as those of max- 
imum density. Even at this stage, we find a qualitative 



resemblance between the Fig ( 10 ), the disk packing prob- 



lem, and the corresponding sphere packings (as shown in 
Fig©). 

A simple method for achieving a semiquantitative de- 
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FIG. 10. Area fraction 11 of disk pacltings on the plane which are consistent with wrapping onto a cy Under of diameter |V|/7r. 
Asymmetric packings are shown by (red) dashed curves; line-slip solutions are (black) continuous curves. Symmetric structures 
labelled [/, m, n] correspond to points of maximum density while intermediate asymmetric packings are labelled using the 
rhombic notation [p, q] 
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FIG. 11. A plot of Eq. Q for various values of D jd, show- 
ing the approach towards elliptical, and eventually circular, 
packing as D jd increases beyond 2. 



scription of the sphere packing problem relies on the fact 
that for the range 1 < D < 2.71486 all the spheres in 
the packing are in contact with the cylindrical boundary. 
Thus there exists an inner cylinder on which the centres 
of all the spheres are to be found. The pattern formed by 



the intersection of the spheres with the surface of the in- 
ner cylinder bears a close resemblance to the disk packing 
problem, which we now consider in detail. 



B, Relation to sphere packing 

We now explore more closely the quantitative corre- 
spondence between disk and sphere packings, and find a 
way of bringing them as close to each other as is possible 
without undue complication. 

Consider again a packing of spheres, all touching the 
confining cylinder of diameter D. Their centres lie on an 
inner cylinder of diameter D' — D — d, so that a given 
centre is located at {D' /2, 6, z) in cylindrical polar coor- 
dinates. The separation of contacting spheres in 3D is d 
but if the inner cylinder is rolled out onto the plane, i.e. 
a given sphere centre is mapped to the 2D cartesian co- 
ordinates (s = D'9/2,z), then the distance between the 
centres of contacting spheres in 2D depends on their mu- 
tual orientation. 

In general, the required contact rule, in 3D, for a pair of 
spheres (each having diameter d) whose centres, in cylin- 
drical polar coordinates, are located at ri = [D' /2, 9, z) 
and ra = {D' /2, 6 + M,z + Az) is [3], 
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D'2 

(Az)2 + (1-C0SA6') 



(4) 
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We plot the parametric equation Eq. Q in Fig 
for various values of D/d (where the distance As 
{D' /2)A9 is referred to as the 'arc length'). We see that, 
as D/d increases above 2, the 'surface packing' resembles 
the packing of oval disks, which may be approximated by 
ellipses. For contacting spheres that lie directly in a line 
parallel to the direction of the cylinder axis, the separa- 
tion of their centres is still 5*1 1 = d, while for those that 
lie in the perpendicular plane, it is, 



S\ ^ D' sin" 



(5) 



For simplicity we take the major and minor axes of the 
ellipse to have the Sji_ and 5*1 1 , so that the distance on the 
plane for a pair of contacting spheres is approximated by. 



S"^ =d^s\u^{(i)) 



D' sin" 



D' 



cos^ ((/)). (6) 



Comparing the circumference of the inner cylinder ttD' 
with the length of the periodicity vector Vd (measured 
in units of the disk diameter d), gives a stretch factor, 



X = 



ttD' 
Yd 



d ' 



(7) 



where the last term is the ratio of the arc distance be- 
tween contacting spheres, which lie perpendicular to the 
direction of the cylinder axis. In this way the sphere 
packing is related to the planar packing of elliptical disks, 
which is merely a stretched version of a circular disk pack- 
ing. Given any packing of circular disks in the plane, it 
may simply be stretched in the direction of the vector V, 
by a factor X, and wrapped onto a cylinder of diameter 



D' 



sin(7r/y) 



(8) 



(obtained using Eq. IT])), to create a good approximation 
to a sphere packing. 

Thus a correspondence that at first seems distant may 



be brought much closer. Fig ( 12 ) presents the resulting 



transformed curves for sphere packing density; asymmet- 
ric and line slip structures are shown by red (dashed) and 
black (continuous) lines, respectively. The peaks in the 
density correspond to symmetric packings. The lower 
(dashed) heavy black line accentuates the curve of max- 
imum density, as predicted by our analytical approach. 
This is to be compared with the results of simulated an- 
nealing shown by the upper curve. We have extended 
our simulated annealing study of single layer packings 
to D > 2.71486, as shown by the grey region in Fig 



(12), by searching for the densest structures in which 



all the spheres have centres on the surface of a cylin- 
der of diameter D' . From this we see that the maxi- 
mum possible density, for single layer packings, is found 



at $(2.4863) = 0.5446 and corresponds to the symmet- 
ric packing (5,3,2). Beyond this the packing density 
steadily diminishes since spheres are to be found only 
on the surface of the cylinder and the interior is empty. 
The analytical method presented here gives the correct 
line-slip solution in the majority of cases. There are, nev- 
ertheless, notable discrepancies between the analytical 
and numerical results. These are due to the simplicity of 
the transformation used to map the disk packing problem 
into the sphere packing problem; as a consequence, the 
analytical results always predict a type (3) line-slip so- 
lutions leading up to n, n, packing while simulations in 
fact find a region split between type (2) line-slip followed 
immediately by a type (3) line-slip. A more accurate (and 
more complicated) transformation ought to account for 
this by pushing the type (2) curve above the type (3) 
solution, for part of this region. This is borne out by the 
fact that as the diameter of the cylinder increases, and 
higher order corrections to the transformation diminish 
in importance, the type (2) solution in the simulations 
are observed over an increasingly smaller range in the 
split region. 



VIII. STRUCTURES BEYOND D/D = 2.71486, 
WITH INTERNAL SPHERES. 



We are able to pursue the same simulations to higher 
cylinder diameters, finding the anticipated occurrence of 
internal spheres, initially at D/d = 2.71486. The struc- 
tures involved are tabulated in the final section of Table 
HI In some cases phyllotactic notation may still be use- 
ful, but much of the simplicity of the previous sections 
is lost, and increasingly so. The extension of the curve 
of volume fraction, as far as present computing resources 
allow, is shown on the right hand side of Fig pi. 

It is of the same general character as before; that is we 
observe peaks corresponding to maximal contact struc- 
tures - which are labelled using the notation (m,n) in 
table I - and decreasing or increasing D/d yields struc- 
tures with fewer contacts which can in some way be de- 
formed steadily into the maximal contact packing - these 
are labelled —{m,n) and +{m,n) respectively. 



Packings with internal spheres are depicted in Fig ( 13 ) 



where in each case the middle diagram shows the en- 
tire structure viewed side-on. Although the introduction 
of inner spheres produces more complicated structures 
there remains an outer shell of spheres contacting the 
confining cylinder. The rolled-out pattern made by these 
outer spheres is shown to the left of each packing. The 
arrangement of the inner spheres is depicted in the di- 
agram to the right of each packing, whereby the outer 
layer is transparent and the inner spheres are in yellow. 
The first packing with internal spheres is structure 30, 
or —(1,5), which consists of a basic unit cell containing 
six spheres. Of these one sphere is not in contact with the 
cylindrical surface and is found at the centre of the cylin- 
der. The remaining five spheres, which are in contact 
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FIG. 12. Comparison of simulated and analytic volume fractions for the densest packings - upper and lower curves, respectively. 
Dotted lines are a guide for the eye to the detailed correspondence. In the analytic results, the line-slip structures are identified 
by black continuous lines, and the red dashed-dotted lines denote asymmetric packing structures. The numbers 1, 2, and 3 (see 
text) denote the type of line slip observed (upper curve) or predicted (lower curve)with 1\2 denoting a degenerate case. The 
shaded region on the right is the continuation of the single layered structures for D/d > 2.71486. 



with the cylinder, form a tilted pentagonal ring around 
the central sphere (i.e. tilted at an angle with respect 
to the plane perpendicular to the cylindrical axis). The 
extended structure can be made by translating the unit 
cell along the cylinder and rotating it by tt. The result is 
a chiral packing with an internal chain of non-contacting 
spheres along the cylindrical axis and an outer layer of 
touching pentagonal rings. As can be seen, in the corre- 



sponding rolled out diagram of the outer layer in Fig ( 13 1 , 
there is only a single point of contact between outer rings 
corresponding to a sphere with four surface neighbours. 
The location of this sphere alternates by an azimuthal 
angle of tt between successive layers. 

Some of the spheres in the outer layer have only two 
contacts with other spheres in the outer layer. However, 
mechanical stability is assured by their contact with the 
internal chain of spheres and the average number of con- 
tacts per sphere for the structure as a whole is above 
four. 

Increasing D jd reduces the tilt of the pentagonal rings 
until at D I d = 2.7306 we have the maximal contact pack- 
ing (1,5). Thus in structure 31 the spheres in the inter- 
nal chain are in contact with each other and, compared 
with structure 30, there are now a greater number of 
contacts between successive pentagonal rings. Increas- 
ing Djd breaks some of these surface contacts and forces 



the chain of internal spheres to form a twisted zigzag 
structure. Thus structure 32, or -|-(1,5), is chiral and is 
superseded at Djd — 2.74804 by a new type of packing. 

Structure 33, or —(1,10) consists of a basic unit cell 
containing eleven spheres. Of these ten touch the cylin- 
drical surface and are arranged as a pair of pentagonal 
rings stacked on top of each other. The eleventh, inter- 
nal, sphere is found above the ten surface spheres and is 
located at the centre of the cylinder. The extended struc- 
ture can be described as an alternating sequence of sur- 
face spheres followed by an internal sphere. However, on 
average each sphere is only in contact with 40/11 = 3.636 
other spheres, this number is too low for mechanical sta- 
bility, which we explain as follows. 

As shown by the rolled out diagram of the outer sphere 
centres, surface spheres from adjacent unit cells are not 
in contact. Thus the pairs of pentagonal rings are free 
to rotate, by a certain amount about the cylindrical axis 
while internal spheres remain fixed in position. Whether 
structure 33 is achiral or chiral depends on the relative 
orientation of the surface spheres with respect to each 
other. No other optimal cylindrical packing yet discov- 
ered has this property. 

With increasing Djd the pentagonal rings from ad- 
jacent unit cells are eventually locked into position at 
Bjd = 2.8211 to give (1, 10). Thus packing 34 is a max- 
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FIG. 13. Structures with internal spheres. In the case of structures which are not maximal contact the images are produced 
from numerical results for arbitrarily chosen values of D/d within their ranges. The central image shows each structure viewed 
side-on. The image on the right shows the inner layer by making the outer layer transparent. The image on the left, for 
each packing, is the rolled out structure generated by the outer layer of spheres which contact the cylindrical boundary - the 
periodicity vector V is indicated by the red arrow. 



iinal contact achiral structure, corresponding to a peak 
in Fig (|3|. The spheres in the outer layer form a perfect 
rhonibic[5, 5] lattice when rolled out onto the plane and 
the internal spheres lie along the cylindrical axis. An 
increase in D/d results in the loss of two contacts (on 
average) between the surface spheres and the internal 
spheres, and +(1, 10) proceeds downwards in density. As 
the trend continues we observe a decrease in the separa- 
tion of neighbouring internal spheres and a modulation in 
the rhombic surface pattern; this latter symmetry break- 
ing is responsible for structure 35 being a chiral structure. 

Packing 36 includes an internal linear chain of non- 
contacting spheres lying along the central axis of the 
confining cylinder. These are surrounded by an outer 
layer of spheres which form a complex chiral structure. 
As D/d is increased the spheres in the internal chain are 



brought ever closer to each other until they make con- 
tact at D/d ~ 2.8615, which corresponds to the maximal 
contact, chiral, structure 37. 

A further increase in D/d forces the internal chain of 
spheres to form a twisted zigzag structure. As a result 
there is a loss of two types of contacts: the breaking of 
contacts between spheres in the outer layer (as seen from 
the corresponding rolled out pattern in Fig (13)) and a 
break in contact between spheres in the outer and inner 
layer. Thus there is a decrease in density as we proceed 
from structure 37 along 38, i.e. as we increase D/d for 
+(2,5). 

Structure 39 is remarkable in that we find an internal 
chain of spheres (along the central axis of the cylinder) 
but the chain is composed of a pair of touching spheres 
followed by a gap, this is the structure —(2, 13). Increas- 
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ing D /d produces structure 40, a maximal contact chiral 
packing (2, 13), which compared with structure 39 has 
an increased number of contacts between internal-surface 
and surface-surface spheres. 

This remarkable sequence of structures would have 
been difficult to imagine in advance. It is likely to be fol- 
lowed by an equally rich scenario, whenever it becomes 
possible to pursue higher values of D/d. The structures 
reported are new, to our knowledge, except in so far as 
they correspond to dry foam structures mentioned below 
in section 10. 



IX. RELATED EXPERIMENTS 

Columnar sphere packings appear in a variety of dif- 
ferent experimental contexts, such as the packing of Cgo 
buckyballs inside carbon nanotubes [TB] and polystyrene 
spheres inside the pores of silicone membranes jl4j . The 
structures that are commonly identified are the straight 
chain (structure 1 in our table I), zigzag (structure 2 
or Ci), twisted zigzag (structure 3) and structure 5, al- 
though more complicated structure have also been ob- 
served, e.g. [O] . 

However, care needs to be taken when directly com- 
paring these experimental observations with the simple 
hard sphere simulations described here, which might well 
only serve as a first guide of what structures to expect. 
The minimization of interaction energy may replace the 
maximization of density as the guiding principle. 

For example the polystyrene particles in the experi- 
ments of 14J have been charge-stabilised and thus repel 
each other. Unlike the situation in our simulations, this 
results in a preference to sit near the wall of the pores 
(i.e. the cylinder wall). 

Also the comparison of our simulation results with 
packings of buckyballs might be limited. Experiments by 
PB] show that while it is possible to fill Cgo into double- 
walled nanotubes down to the ratio D/d =1.13 (resulting 
in a simple linear chain), it is not possible to fill single- 
walled carbon nanotubes below D/d = 1.25 {D is the 
inner diameter of the tube in both cases) . This reflects 
the role of the van der Waals forces between Ceo and the 
confining nanotube walls. 



Present experiments with ball bearings and 
bubbles 



However, we have identified two other experimental 
systems that provide very good comparison with our sim- 
ulation data (details to be provided in a follow-up paper) . 

The first are metal spheres of a few millimeter in di- 
ameter (the type used in "ball bearings") packed into 
perspex tubes. The filled tubes were mechanically agi- 
tated over an extended period. At the end of this "an- 
nealing" procedure this resulted in clearly identifiable or- 
dered sphere packings. These experiments were carried 



out for a large range of values of D/d (up to 3.15), and 
resulted in about 16 different structures, of both the max- 
imal contact and the line-slip type. 

We have also carried out more extensive experiments 
with equal-volume gas bubbles (with diameter of a few 
hundred microns). These are produced by bubbling gas 
into surfactant solution with the bubble diameter easily 
variable by adjusting the gas flow. 

The bubbles are then gathered into vertically placed 
cylinders. They maintain a (nearly) spherical shape, even 
when in contact, up to a column height of a few millime- 
tres, corresponding to the capillary length for the surfac- 
tant solution in use. 

We were able to produce straight chain, zigzag and 
a large number of maximal contact structures, and also 
many structures with internal bubbles, up to and well be- 
yond the point reached by the present simulations. There 
is also some evidence of the twisted zigzag structure and 
line-slip structures. The subsequent paper presenting 
these experimentally obtained structures will include the 
use of X-ray tomography to establish their internal con- 
figurations. 



B. Ordered dry foam structures 

Extensive and detailed results have been published for 
the various structures that form when equal volume bub- 
bles with diameter exceeding a few millimeters (i.e. larger 
than the capillary length) are collected in vertical tubes 
[HI [32] . In this case, the liquid drains and a dry foam is 
created - a packing of polyhedral bubbles, with its volume 
fraction $ approaching unity. 

The surface pattern displayed by these ordered foams 
is hexagonal (apart from the so-called bamboo structure 
which is simply an array of parallel soap films). The 
phyllotactic notation is thus the obvious choice for their 
classification, at least as long as there are no internal 
bubbles (see for example Fig. 14). 



Restricting ourselves for the moment to this case it 
would appear that all ordered dry foam structures have 
corresponding sphere packings (wet foams) , i.e. they can 
be classified by the same phyllotactic notation. These 
include the straight chain (bamboo), and zigzag struc- 
ture (called 2-1-1 or staircase structure in the foam lit- 
erature), and all of the 10 maximal contact structures, 
except structure 25 (5-4-1) which has not yet been ob- 
served. None of the ordered dry foam structures is of the 
line-slip type. 

Note that since the bubbles are deformable, the aver- 
age contact number in the case of dry foams is generally 
higher than in hard sphere packings. 

Unlike their sphere packing counterparts, the various 
dry foam structures are found over ranges of values of 
D/d (in the foams literature d refers to the equivalent 
sphere diameter of a bubble). Hysteresis plays a large 
role in the standard experimental procedure, leaving to 
overlapping ranges of stability for the various structures. 
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(a) 



(b) 



FIG. 14. In the dry foam equivalent of the zigzag structure 
(often called the staircase structure) equal-volume gas bubbles 
are separated by thin liquid films, (a) Computer simulation 
using the Surface Evolver Software of Ken Brakke 23^. (b) 
Photograph of bubbles in a cylinder of about 1cm in diameter. 



The structure of minimal energy for a given value of 
D /d may be determined from computer simulations [221 
f24U27| . We find these ranges always to be lower than 
the value of D/d for the corresponding sphere packing, 
reflecting their much smaller volume fractions. It remains 
for future research to establish the phase diagram (which 
may be rather complex) that connects these two limiting 
cases. 

There are also observations of over 20 different ordered 
dry foam structures with internal bubbles [8j , but the pre- 
cise internal arrangements have only been identified for 
the three simplest cases, with the aid of Surface Evolver 
simulations |23J[2S]- 

The dry foam with a surface pattern equivalent to 
structure 28 (C5) has six bubbles in the periodic unit 
cell, i.e. one more than C5. A pentagonal dodecahedron 
in the centre of the cylinder is surrounded by a ring of 
five bubbles in contact with the cylinder surface. The 
dodecahedra of neighbouring unit cells are in contact. 

There is also a foam structure with a so-called Kelvin 
cell (tetracaidecahedron) in the centre, surrounded by 6 
bubbles in contact with the surface. Again the internal 
bubbles of neighbouring unit cells are in contact with 
each other. 

The third identified structure with internal bubbles 
consists of a a total of 13 bubbles in the unit cell, 12 
touching the surface, and a bubble with 12 pentagons 
and 3 hexagons (Goldberg-3) in the centre [5S]. Inter- 
estingly in this case the internal bubbles of neighbouring 
unit cells are not in contact, similar to the also achiral 
structure 34. The surface pattern of this foam structure 
consists of an arrangement of bi-disperse hexagons. 



X. CONCLUSIONS 

Our simulations have provided detailed results for 40 
distinct columnar crystals, of which many are new. They 
fall into three categories: the very simplest cases, where 
the sphere diameter is of the same order as that of the 
cylinder, a wide range of further "phyllotactic" struc- 
tures, which may be understood as related to surface 
packings of disks, and structures that incorporate inter- 
nal spheres. 

Corresponding experimental observations are now 
available and will be reported in a second paper. Ex- 
periments can already be pushed to much larger cylinder 
diameters, providing further insights and challenges to 
simulation and interpretation. 

Such a further development should motivate a reassess- 
ment of the simulation methods to be used; we make no 
strong claims for the efficiency of that used here. We have 
necessarily been cautious in using it, stipulating high 
degrees of convergence and undertaking many repeated 
runs in search of the optimum, to increase confidence in 
our conclusions. In parallel studies, Chan [3] used a se- 
quential deposition method that was highly expeditious 
in determining optimal structures within the restriction 
of that procedure, whose effect could not be known a 
priori. Suitably adapted, it was able to reproduce the 
structures reported here, up to at least D/d — 2.7013. 
Future work will include an extension of this deposition 
approach to higher values of D/d. 

Finally we note some possible extensions to this work. 
An example is the study of disordered hard sphere pack- 
ings in a cylinder; for such problems the definition of the 
volume fraction recently given by Chan [3] may prove use- 
ful. A simple problem that may be of significant practical 
interest is to find the densest packing of hard spheres in 
a cylinder which is capped on either side by hard walls 
(i.e. packing in a cylindrical box). Even richer possi- 
bilities are offered by packing hard spheres in a cylinder 
that is capped on both ends by spheres held in a fixed 
position. Indeed, hard sphere packings that are bounded 
by templated surfaces on all sides are of significant cur- 
rent interest and are capable of realising a range of dif- 
ferent morphologies (including the Weaire-Phelan struc- 
ture) [IH]- 
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XII. APPENDIX A: SIMULATION TECHNIQUE 
A. Energy function 

The simulation is addressed to a cylindrically shaped 
cell of length L and diameter D. Contained within this 
space are N points which represent the centres of N 
spheres, each of diameter d. If a pair of spheres is suffi- 
ciently close that they overlap we account for this using 
a pairwise potential as described below. A similar over- 
lap potential is used to prevent the spheres from escap- 
ing the simulation cell in the radial direction. The final 
structures will be the densest that we find which have 
zero overlap energy. 

In addition, we impose twisted periodic boundary con- 
ditions on the top and bottom of the cylindrical cell as 
described in the main text, with a twist angle a. 

We model the overlap potential between spheres us- 
ing a Hookean, or "spring-like" , pairwise interaction be- 
tween the ith and jih. spheres, which have their centres 
at r^ = {ri,9i,Zi) and Vj — {rj,9j, Zj), the interaction 
energy between spheres is then given by. 



Er, 




df 



if Ti 
if r,; 



<d 
>d 



(9) 



where r, 



r,| is the distance between the centres 



of the spheres. Note the interaction energy falls to zero 
when there is no overlap between the spheres. 

The interaction energy between the ith sphere and the 
boundary is given by. 



Ef 



UnB-d/2y 



if r,B < d/2 
if ne > d/2 



(10) 



a by using a random number generator. A small initial 
value for the cylinder length L is chosen to insure overlap 
between the spheres. 

Keeping D and L fixed, we search for the lowest energy 
arrangement for the N spheres by varying their coordi- 
nates and the twist angle. This is done using the standard 
Metropolis simulated annealing algorithm, where for a 
cluster of N spheres the algorithm was run with typically 
iV X (5 X 10^) Monte Carlos steps. The temperature of the 
simulation was decreased linearly. The average displace- 
ment of the spheres at each temperature step was chosen 
by an automatic process to give an acceptance probabil- 
ity of 0.5 ± 0.01. The results of the simulated annealing 
are then put through a conjugate gradient routine to en- 
sure that a local minimum has been reached. 

The whole process is repeated many times, using a new 
randomly generated initial configuration, to give confi- 
dence that the lowest energy state has been found. From 
this ensemble we take the lowest energy configuration as 
the final state for that particular run. 

After this first run we perform a subsequent run with 
a slightly longer cylinder. Since the spheres have more 
room the energy of the final state is lower compared to 
the final state of the previous run. Using a divide and 
conquer approach we are able to establish the cell length 
at which the energy per sphere in the system falls to zero, 
in practice this means the value of the energy is within a 
critical bound which is set to be E/N=l x 10^^±5 x 10^^. 
At this point the spheres have a small overlap correspond- 
ing to a five decimal place accuracy in the volume fraction 
(this is deduced by comparing simulation with the analyt- 
ically derived volume fraction for the Cn circle packing 
structures). 

From this final result we compute the volume fraction 
which is defined as, 



where r^s = \D/2 — ri\. 

The twisted periodic boundary conditions on the ends 
of the cylinder are incorporated as follows: The ith 
particle in the simulation cell has an image at the top 
and bottom of the simulation cell, the coordinates of 
these images are given by rf — {ri, 9i -\- a^Zi + L) and 
r~ = (r^, 9i — a, Zi — L), respectively, where L is the length 
of the cylinder and a is the twist angle. 

Thus the total energy of the system is given by the sum 
of the sphere-sphere, sphere-boundary and sphere-image 
interactions. 



^{N,D) 






(11) 



where Vg = {A/3)TT{d/2)^ is the volume of one of the 
hard spheres and Vc — tt{D/2)'^L is the volume of the 
simulation cell. Note that the volume fraction depends 
on both D and N since different values of N yield unit 
cells with different structures. Thus the whole procedure 
is repeated for series of different values of A^ = 1, 2, 3..., 
typically up A'^ = 15 for the structures without internal 
spheres, before accepting the one which gives the largest 
volume fraction, as shown in Table IT] 



B. Numerical Method 

For a tube of diameter D we wish to find the unit 
cell, composed of A^ spheres, which when rotated and 
stacked along the tube has the highest volume fraction. 
In this section we describe the simulation protocol used 
to achieve this. 

For a given D and d we assign initial starting positions 
to the A^ spheres and an initial value to the twist angle 



XIII. APPENDIX B: CONTINUITY OF THE 
VOLUME FRACTION 

We demonstrate that there can be no sudden finite dis- 
continuities in the curve describing the maximum volume 
fraction (density) as D/d is varied. This follows from 
lower and upper bounds, as described below, which limit 
the variation in volume fraction in the neighbourhood of 
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any point Dg/do. We will see that square root singu- 
larities observed in the numerical results for the volume 
fraction are due to the arguments given below. 



A. Low^er bound 

For increasing D/d an obvious variational argument 
bounds the density below. As D/d is increased from 
Dp/do, the cylinder expands radially. We may take the 
structure which has maximum density at Dg/do as a trial 
structure for D/d > Do/do. The volume fraction of the 
structure, at D^/dQ is, 



HDo/do) = 



V, 



^(^o/2)2L' 



and that of the trial structure is, 



^^{D/d) 



V, 



7r(D/2)2Z' 



(12) 



(13) 



where Vg = Vs (d) is the volume of a sphere of diameter d 
and L is the average separation between spheres in the z 
direction. Since $'^(Do/do) = ^{Do/do) it follows that, 



<^^{D/d) > ^^Do/do) 
providing a lower bound for D > Do- 



(14) 



structure has the volume fraction. 



^TiD/2)^L + A) 



(16) 



which when combined with Eq. (12) results in the fol- 
lowing bound 



^'^{D/d)<^j^^Do/do 



(17) 



Hence we again arrive at a (lower) bound for ^^{D/d) 
which goes continuously to ^{Do/do) as D/d — ?► D^/do 
but in this case with a square root form. 



B. Upper bound 



For decreasing D/d, a more subtle argument bounds 
the variation of density below by a square-root function. 
Decreasing D/d forces the spheres to move radially in- 
wards to avoid contact with the cylindrical boundary. 
The resulting overlap between spheres can be eliminated 
by displacement of their centres parallel to the cylinder 
axis. 

Let us index the sphere centres in ascending height us- 
ing the index j - so that Zj > Zj^i for j > j — 1. Given 
the structure at Do/do which has a maximum density, 
then an extreme case is one where successive spheres all 
have the same height. Consider in this case a pair of con- 
tacting spheres with separation d. Reducing the cylinder 
diameter by a factor of X, so that the new diameter is 
D = DqX, will force an overlap so that the separation 
between sphere centres is now dX. The overlap can be re- 
moved by moving one of the spheres vertically a distance 
A, so that their separation is once again d. It follows 
that, d"^ = A^ — (dX)^. Thus a constant C can be chosen 
so that the following choice is sufficient. 



A = CjD^^ 



D^ 



(15) 



In order to eliminate the overlap, the sphere centres are 



shifted to a new height z' 



Zj + 



jA. The resulting trial 
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